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Abstract 

A class of A.L.E. time discretisations which inherit key energetic properties (non- 
linear dissipation in the absence of forcing and long-term stability under conditions 
of time dependent loading), irrespective of the time increment employed, is estab- 
lished in this work. These properties are intrinsic to real flows and the conventional 
Navier-Stokes equations. 

A description of an incompressible, Newtonian fluid, which reconciles the differ- 
ences between the various schools of A.L.E. thought in the literature is derived for 
the purposes of this investigation. The issue of whether these equations automat- 
ically inherit the afore mentioned energetic properties must first be resolved. In 
this way natural notions of nonlinear, exponential-type dissipation in the absence 
of forcing and long-term stability under conditions of time dependent loading are 
also formulated. 

The findings of this analysis have profound consequences for the use of certain 
classes of finite difference schemes in the context of deforming references. It is 
significant that many algorithms presently in use do not automatically inherit the 
fundamental qualitative features of the dynamics. 

The main conclusions are drawn on in the simulation of a driven cavity flow, a 
driven cavity flow with various, included rigid bodies, a die-swell problem, and a 
Stokes second order wave. The improved, second order accuracy of a new scheme 
for the linearised approximation of the convective term is proved for the purposes of 
these simulations. A somewhat novel method to generate finite element meshes au- 
tomatically about included rigid bodies, and which involves finite element mappings, 
is also described. 



Keywords: Energy conservation; incompressible, newtonian fluid; completely general ref- 
erence description; arbitrary Lagrangian Eulerian; A.L.E.; free surface; finite elements; 
new Poincare inequality; second order accurate linearisation of the convective term; au- 
tomatic mesh generation. 
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1 Introduction 

This work focusses on establishing a class of A.L.E. time discretisations which inherit key 
energetic properties (nonlinear, exponential-type dissipation in the absence of forcing and 
long-term stability under conditions of time dependent loading) irrespective of the time 
increment employed. The findings of this analysis have profound consequences for the 
use of certain classes of difference schemes in the context of deforming references. It 
is significant that many algorithms presently in use do not automatically inherit the 
fundamental qualitative features of the dynamics. 

Descriptions of fluid motion are conventionally based on the principles of conservation of 
mass and linear momentum. One might hope that all such descriptions would accordingly 
exhibit the afore mentioned, key energetic properties consistant with the principle of 
energy conservation. These properties are intrinsic to real flows and the conventional, 
Eulerian Navier-Stokes equations. 

A description of an incompressible, Newtonian fluid, which reconciles the differences be- 
tween the so-called arbitrary Lagrangian Eulerian (A.L.E.) formulation of Hughes, Liu 
and Zimmerman [|] (deformation gradients absent) and that of Soulaimani, Fortin, 
Dhatt and Ouellet [|r3| (deformation gradients present, but use is problematic), is de- 
rived for the purposes of this investigation. The implications of the resulting description 
are investigated in the context of energy conservation in a similar, but broader, approach 
to that taken by others (eg. SiMO and Armero [|I4|) for the conventional, Eulerian 
Navier-Stokes equations. 

The main conclusions of this work rely on a new inequality and a number of lemmas, the 
proofs of which are listed in an appendix at the end of the paper. The new inequality is 
used in place of where the Poincare-Friedrichs inequality might otherwise have limited the 
analysis. The lemmas are mainly concerned with the new convective term. This analysis 
is extended in that non-zero boundaries, so-called free boundaries and time-dependent 
loads are considered. 

The resulting theory is used in the simulation of a driven cavity flow, a driven cavity flow 
with various, included rigid bodies, a die-swell problem, and a Stokes second order wave. 
A new scheme for the linearised approximation of the convective term is proposed and 
the improved, second order accuracy of this scheme is proved for the purposes of these 
simulations. A somewhat novel method to generate flnite element meshes automatically 
about included rigid bodies, and which involves finite element mappings, is also described. 



2 A Completely General Reference 

The implementation of most numerical time integration schemes would be problematic 
were a conventional EulerianQ description of fluid motion to be used in instances involving 
deforming domains. The reason is that most numerical time integration schemes require 



Eulerian or spatial descriptions are in terms of fields defined over the current configuration. 
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successive function evaluation at fixed spatial locations (the exception being the finite 



element with respect to time approach of Tezduyar, Behr and Liou [|T8[). On the 
other hand meshes rapidly snarl when purely LagrangianQ descriptions are used. It is for 
these reasons that a completely general reference description is usually resorted to. 

Eulerian and Lagrangian references are just two, specific examples of an unlimited number 
of configurations over which to define fields used to describe the dynamics of deforming 
continua. They are both special cases of a more general reference description, a descrip- 
tion in which the referential configuration is deformed at will and which is the focus of 
this investigation. A deforming finite element mesh would be a good example of just such 
a deforming reference in practice. The transformation to the completely general reference 
involves coordinates where used as spatial variables only and the resultant description is 
therefore inertial in the same way as Lagrangian descriptions are. 



2.1 Domains, Mappings and a Notation 

Consider a material body which occupies a domain Q at time t. The material domain, 
Qq, is that corresponding to time t = to (the reference time, to, is conventionally, but 
not always, zero). A third configuration, Q, which is chosen arbitrarily is also defined for 
the purposes of this work. The three domains are related in the sense that points in one 
domain may be obtained as one-to-one invertible maps from points in another. 

For any general function f{x,t), a function, f{x,t) = f{X*{x,t),t), can be defined in 
terms of the domains and one-to-one, invertible mappings illustrated in Fig. ||. Similarly, 
fo{xo,t) = f{X{xo,t),t) can be defined. This notation can be generalised for the 
component- wise definition of higher order tensors. The key to understanding much of 
this work lies possibly in adopting a component-wise defined notation. 

In contrast to the function notation just established, the definition of the operators V 
and div is not based on V and div. They are instead the referential counterparts, that is 

~ d . d d d 

V = 7— and div = 7— — h 7— — h 



dx dxi dx2 dx 



3 



The notation A : B is used to denote the matrix inner product AijBij throughout this 
work, ( ■ , ■ )/,2( . ) denotes the inner product and || ■ ||i2( . ) the norm. 



2.2 Some General Results for Functions Defined on the Three 
Domains 

Three important results are necessary for the derivation of the completely general refer- 
ence description and these are presented below. 



^Lagrangian or material descriptions are made in terms of fields defined over a reference (a 
material reference) configuration. 
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•X ] 

^ ] 

y Spatial Domain 

.-^ (in terms of which the 

conventional Eulerian 
description is made) 




{Si is spatially distorting Si, a 
deforming mesh, for example.) 

Figure 1: Schematic Diagram of Domains and Mappings Used in a completely general 
reference Description 

The Material Derivative in Terms of a Completely General Reference 

The material derivative of any vector field v in terms of a completely general, reference 
is 

(1) 

where is the velocity of the reference deformation, and F is the deformation gradient 
given by 



F{x) = 




This result is demonstrated in Appendix II. 



Material Domain 




dv 
'dt 



F {v 
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An Element of Area in Terms of a Distorting Reference 



The second important result can be recalled from general continuum mechanics. Consider 
an element of area, size dA, with an outward unit normal n. Then 

ndA = F'^NJdA (2) 

where dA and N denote the respective analogous size and outward unit normal of this el- 
ement of area in the referential configuration and J = det F. This result is demonstrated 
in most popular textbooks on continuum mechanics (eg. Lai, Rubin and Krempl ||10|| ). 

Remark: Notice that N is the single exception to the ~ quantities devised in this work. 
A?" is the surface normal perceived in a completely general reference and the components 
of n and N need have nothing in common [N is not n). 



The Kinematic Relation J'q = J'q div v 

The material derivative of the Jacobian j7o is given by the relation 

Jo = Jo div V 

where Jo is defined as follows, 

Jo = det 



dxo J 

This result is demonstrated in most popular textbooks on continuum mechanics (e^ 
Marsden and Hughes |n[]). 



2.3 Derivation of the Completely General Equation 

One way in which to derive a completely general reference description of an incompress- 
ible, Newtonian fluid is to start with the balance laws in global (integral) form, and 
to make the necessary substitutions in these integrals. The desired numerical imple- 
mentation (similar to the conventional Navier-Stokes one which has been thoroughly 
investigated and found to be stable) is then obtained. 



Conservation of Mass 

Let Q{t) be an arbitrary sub- volume of material. The principle of conservation of mass 
states that 

d f 

— / pdQ = (rate of change of mass with time = 0) 
dt Jn{t) 

d f 

-r / PoJod^o = (reformulating in terms of the material 
at JUn 
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configuration, VLq.) 

d 



f o 

/ TT {Po>^o} dflo = (since limits are not time dependent in 
Jno ot 

the material configuration.) (3) 

J {pojid + PqJo) dVlo = (by the chain rule) 

f ip + p div v)dn = (us^ng the k^emat^c result J, = Jodiv v) 
Jn(t) 



in{t) 

dvi dxj 



/ P + Ptt^tt^ Jd^ = (reformulating in terms of the distorting 
Jn{t) \ oxj oxi ) 



referential configuration, Vt{t).) 
p + pSJv : F J = (integrand must be zero since the volume 

was arbitrary.) 

Thus, for a material of constant, non-zero density, 

Vv : F =0 since J 7^ (mappings are one-to-one and invertible) . 

Notice also that equation @ implies 

^{PoJo} = (4) 
since the volume was arbitrary and the integrand must therefore be zero. 

Conservation of Linear Momentum (and Mass) 

The principle of conservation of linear momentum for an arbitrary volume of material 
fl{t) with boundary T(t) states that 

— / pvd^l = [ phdVt + / (TndA (5) 

at Jn{t) Jn{t) Jr(t) 

where p is density, b is the body force per unit mass, cr is the stress, n the outward unit 
normal to the boundary and v is the velocity. The term on the lefthand side can be 
rewritten as follows: 



— / pvdQ = — poVaJ'^dflo (Reformulating in terms of the material 
at Jn(t) at Jno 

configuration, Qq.) 

f d 

— ^ {PqVoJ^o} dQo (Since limits are not time dependent in 

Jflo ot 

the material configuration.) 
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Q(t) 



pvdVt 



pvJdVt 



n(t) 



f dv 

Q{t) ^ \dt 



(The second term above is zero as a 

consequence of equation 
(Reformulating in terms of the dist- 
orting referential configuration, Q.) 



Vv F 7v — v"'"-') ] JdQ. (Using result 

(0; on page ^ 

where i) denotes the material derivative of v. The surface integral becomes 



crndA 



(Reformulating in terms of a distorting 
reference using result on page \^.) 



div {cri^ J^dVt (By the divergence theorem). 



n{t) 

Finally, the term involving body force becomes 



phdVt 



pbJdCl 



(Reformulating in terms of a distorting 
reference.). 



Substituting these expressions into (|^), remembering that the volume used in the argu- 
ment was arbitrary and that the entire integrand must therefore be zero, the conservation 
principles of linear momentum and mass may be written in primitive form as 



dv 



pi^ + VvF \i-v'''f)\j = pbJ + divP 



and 



Vi) : F * 



(6) 



(7) 



where P is the Piola-Kirchoff stress tensor of the first kind, P = a-F J. In terms of 
the constitutive relation, cr = —pi + 2/iZ), for a Newtonian fluid. 



-pi + p 



VvF + VvF 



F-'j 



smce 



D = - [Vv + [Vv 



The derivation of a variational formulation is along similar lines as that for the Navier- 
Stokes equations (the purely Eulerian description). For a fluid of constant density, the 
variational formulation 



p I w ■ ^—JdVt + p I w ■ Vv 



dt 



[V — V 



JdQ 



p I w ■ hJdVl 



pVw : F ^JdVt - 2p I D{w) : J:>{v)Jdn 



wPNdT 
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L 

Jn 



qVv : F *dn 







(9) 



is obtained, where q and w are respectively the arbitrary pressure and velocity of the 
variational formulation. 

Notice that the usual procedure of assigning a value of zero to the arbitrary velocity, 
w, at the boundary has not been followed. The boundary integral in the variational 
momentum equation has consequently not been eliminated as is normally done. The 
reasons are twofold; firstly problems for which the ensuing investigation is intended are 
of a free boundary type and so the solution is not known there; secondly, a specific 
function (which cannot arbitrarily be assigned a value of zero at the boundary) will be 
substituted for w in the forthcoming analysis. 

2.4 Reconciling the Different Schools of Thought 

The equations (|^) and (0) are the completely general referential description of an incom- 
pressible, Newtonian fluid. They reduce to the so-called A.L.E. equations of Hughes, 
Liu and Zimmerman for an instant in which spatial and referential configurations 
coincide. 

Since the approximate set of equations is broken into a sequence of discrete time steps in 
the implementation, one is entitled to choose a new referential configuration during each 
time step, should one so desire. This is what is known as an "updated" approach; when 
each time step is really a fresh implementation. In the case of time stepping schemes 
based about a single instant (eg. the generalised class of Euler difference schemes to 
be investigated in Section ^ a considerably simplified implementation can be achieved 
by an appropriate choice of configurations. Making the choice of a referential configura- 
tion which coincides with the spatial configuration at the instant about which the time 
stepping scheme is based allows the deformation gradient to be omitted from the approx- 
imation altogether (the deformation gradient is identity under such circumstances). For 
such implementations (those which require evaluation about a single point only) no error 
arises from the use of the equations cited in Hughes, Liu and Zimmerman [^, 



These equations are not valid for any, arbitrary choice of reference or if the implemen- 
tation requires the equation to be evaluated at more than one point within each time 
step (eg. a Runge-Kutta or finite-element-in-time scheme). It is important to remember 
that in a discrete context the reference configuration is fixed for the duration of the entire 
time increment. Although the referential configuration is hypothetical and can be chosen 
arbitrarily for each time step, once chosen it is static for the duration of the entire time 
step. Once the coincidence of configurations is ordained at a given instant, F is defined 
by the reference (mesh) deformation, both before and after, and must be consistant. 




pb + div cr 



(10) 



div V 



0. 




The equations of Hughes et al. are an arbitrary Lagrangian Eulerian (A.L.E.) de- 
scription in the very true sense under the circumstances of implementations requiring 
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evaluation about more than one point within each time step (this is not surprising con- 
sidering the equations have their origins in the arbitrarily, either Lagrangian or Eulerian 
programmes of HiRT, Amsden and CoOK 0). This fact is further bourne out in observ- 
ing that key energetic properties, consistant with the principle of energy conservation, 
are not automatically inherited by the equations of Hughes et. al. in the context of 
more general references. 



The momentum equations of SOULAIMANI, Fortin, Dhatt and Ouellet |[TH| are 



flawed as a result of the mistaken belief that aF J is the Piola-Kirchoff stress tensor 
of the first kind (pg. 268 of Soulaimani et al.). Yet another problem is illustrated 
by rewriting the conventional incompressibility condition using the chain rule. The new 
incompressibility condition which arises is most certainly 

dvi dotj dvi dxi 

^TT^ = and not i^tt- = 0- 

OXj OXi OXj OXj 

Further errors arising (eg. J omitted in the first term on the right hand side of the 
momentum equation, equation (10) on pg. 268 of Soulaimani et al.) make the use of 
these equations problematic. 

There would seem to be no reason why one would wish to define the deformation about 
a configuration other than that at the instant about which the implementation is based 
(assuming the implementation used is indeed based about a single point eg. a finite 
difference) thereby involving deformation gradients. Resolving the resulting difficulties 
associated with the deformation gradients by means of a perturbation seems unnecessarily 
complicated in the light of the above reasoning. 



Natural Notions of Energy Conservation in Terms 
of the Completely General Equations 



The effect of quantities parameterising reference deformation on key energetic properties 
- nonlinear, exponential-type dissipation in the absence of forcing and long-term stabil- 
ity under conditions of time dependent loading - is investigated in this section. These 
properties, of form 

lim sup K{v 



K{v) < K{v 



to 



-2yCt 



and 



< 



t;||^2(Q) is the total kinetic energy), are intrinsic to real 
Eulerian Navier-Stokes equations (see Temam |jl6[, ||17|| , 



respectively (where X = |p| 
flows and the conventional, 
CONSTANTIN and FoiAS and SiMO and Armero fl^ in this regard). The effect of 
v'"^'^ on the afore mentioned aspects of conservation of the quantity 



k{v) ^ -p 



L2{Q) 



is essentially what is being investigated, with a view to establishing a set of conditions 
under which the discrete approximation can reasonably be expected to inherit these self- 
same energetic properties. 
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One might anticipate key energetic properties to be manifest only in instances involving a 
fixed contributing mass of material, whether its boundaries be dynamic, or not. An anal- 
ysis of this nature only makes sense in the context of a constant volume of incompressible 
fluid. 

Inequalities of the Poincare-Friedrichs type are a key feature of any stability analysis of 
this nature. Gradient containing terms need to be re-expressed in terms of energy. 
In the case of a "no slip" (f = 0) condition on the entire boundary the situation is 
straightforward, in that it is possible to use the standard Poincare-Friedrichs inequality: 
there exists a constant Ci > such that 



The use of the classical Poincare-Friedrichs inequality is otherwise identified as a ma- 
jor limitation, even in the conventional Navier-Stokes related analyses. The Poincare- 
Friedrichs inequality is only applicable in very limited instances where the value for the 
entire boundary is stipulated to be identically zero. For boundary conditions of a more 
general nature, such as those encountered in this study, in which parts of the boundary 
may be either a free surface, have an imposed velocity or be subject to traction conditions, 
a more suitable inequality is required. (It should, however, be noted that subtracting 
a boundary velocity and analysing the resulting equation is nonetheless still a feasable 
alternative, despite the fact that the equations are nonlinear. This approach requires a 
more sophisticated and involved level of mathematics^.) 

Some common boundary types and associated descriptions are briefly summarised as 
follows: 

1. Fixed impermeable boundaries: The description at such boundaries is usually 
Eulerian and the quantities F and v"^^^ consequently become identity and zero 
respectively. These are typically (but not always) "no slip" boundaries, implying 
that V |r= 0. 

2. Free boundaries: Conventional use allows the spatial mesh to slide along free 
boundaries while still maintaining their overall Lagrangian character. Stated more 



The total volume is nonetheless still a material volume overall. 

3. Imposed velocity-type boundaries: Conventional use entails descriptions which 
usually become pure Eulerian at such boundaries. The total flow across such bound- 
aries is zero for an incompressible fluid if volume is to be preserved. For boundary- 
driven flows one therefore usually assumes that the quantity 



vanishes (a la Temam [|T7| ). 

^Temam succeeds in arriving at an estimate which proves the existance of a maximal attractor 
in two dimensions in this manner. 



v\\l^ < Ci\\Vv\\l2 for all v e [//o(fi)]". 



formally. 
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4. Imposed traction-type boundaries: A variety of descriptions are used at such 
boundaries, ranging from pure Eulerian to tlie vanisliing h- (v — v^*^^^ type described 
for free boundaries. 



These are the modes of reference deformation commonly used at boundaries encountered 
in practice and which will need to be accommodated if the theory is to be applicable. 

The particular types of geometry considered are those that arise in problems involving 
the motion of rigid bodies such as pebbles on the sea bed; thus a free surface is present, 
and the domain may be multiply connected. The Poincare-Friedrichs inequality does, 
furthermore, not hold on subdomains of the domain in question and the constant is not 
optimal. 

Further investigation (communication [Q) reveals a similar relation, the so-called 
Poincare-Morrey inequality, holds providing the function attains a value of zero some- 
where on the boundary. The proof of the Poincare-Morrey inequality is, however, similar 
to that of one of Korn's inequalities (see, for example, KiKUCHi and Oden |^). In 
particluar, it is non-constructive, by contradiction and the constant cannot therefore be 
determined as part of the proof. Viewed in this light the forthcoming inequality amounts 
to a specification of the hypothetical constant in the Poincare-Morrey inequality for 
domains having a star-shaped geometry. 



Inequality 1 (A New "Poincare" Inequality) Suppose v is continuous and dif- 
ferentiable to first order and that v attains a maximum absolute value, c, on an included, 
finite neighbourhood of minimum radius Rmm about a point aj^^s"^ (as depicted in Fig. 
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If Q is a bounded, star-shaped (about a point a:;°"§™jQ domain in R^, then 



3 -Rmax Rmin 

where i?max is the distance to the farthest point in Q from as^'^in^ (Proof in Appendix II.) 

Remark: Notice that "no slip" boundaries which contravene the star-shaped require- 
ment are not of any consequence. This is since additional contributions to inequality 
terms, arising due to the inclusion of any such intruding domains, can be arbitrarily 
costructed to have a value of zero (without any loss of generality). If for example, one 
were to apply the inequality in an investigation of the flow around an aeroplane wing, 
one might make the convenient choice of the wing interior as the desired neighbourhood. 

Remark: Notice that "no slip" boundaries which contravene the star-shaped require- 
ment are not of any consequence. This is since additional contributions to inequality 
terms, arising due to the inclusion of any such intruding domains, can be arbitrarily 
costructed to have a value of zero (without any loss of generality). If for example, one 
were to apply the inequality in an investigation of the flow around an aeroplane wing, 
one might make the convenient choice of the wing interior as the desired neighbourhood. 

This inequality is similar to the Poincare-Friedrichs inequality when c = 0, but is ex- 
tended to a geometrical subclass of domains which have free and partly non-zero bound- 
aries in addition to being more applicable to more challenging examples such as the flow 
around an aeroplane wing. It has a further advantage in that the coefficient is an order 
of magnitude more optimal when used under the "no slip" Poincare-Friedrichs condition 
(under such conditions the domain can always be deconstructed into a number of sub- 
domains in which Rmm = |-Rmax)- The Poincare-Friedrichs inequality is a special case of 
the above inequality. The necessary lemma (below) follows naturally from Inequality |I]. 

Lemma 1 (Deviatoric Stress Term Energy) The kinetic energy satisfies the bound 

C ~ ^ ~ ^ ~i 2 

— K{v) < D(v)j2 ^ , J where C is related to the constant in Inequality C > 0. 
(Proof in Appendix II.) 



As the reviewers rightly point out, this is Korn's inequality with a specified constant 
limited to star-shaped geometries and a variety of such inequahties can be found on page 



323 of Marsden and Hughes [11 |. The following lemma will facihtate the ehmination 



of the convective energy rate in the forthcoming analysis. 
Lemma 2 (Convective Energy Rate) The relation 

-p(e,(ve)F-'(e-*"')j)^^^_^ = -\p{v,v{f-'n ■ {v- v'")) j , ^^^^^ 



1 



-p ( V, V—— ) 
2' \ ' dt I 



L2{Q) 



^by which is meant that every point in the domain can be reached by a straight hne from a;°'^'S™ that 
does not pass outside of Vl 
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holds under the conditions required for equations (jl^ and 
reference description. (Proof in Appendix II.) 

The above lemma is crucial to the analysis for deforming references in particular. The 
following lemma will establish that the boundary term vanishes at free boundaries under 
conditions of conventional usage. 

Lemma 3 (Free Boundary Energy Rate) The boundary term 

vanishes at free boundaries provided the description there is of a vanishing h- (^v — v^^^^ 
type. (Proof in Appendix II.) 

This concludes the preliminaries required for the deforming reference energy analysis. 



11) to be a completely general 



3.1 Exponential Dissipation in the Absence of Forcing 



The issue of whether nonlinear, exponential-type dissipation in the absence of forcing is 
a property intrinsic to the deforming reference description is resolved by the following 
theorem. 



Theorem 1 (Exponential Dissipation in the Absence of Forcing) A sufficient 
condition for the completely general reference description to inherit nonlinear, exponential 
type energy dissipation 

in the absence of forcing is that the reference moves in a vanishing fi- {v — v^^^^ fashion 
at free boundaries and becomes pure Eulerian at boundaries of a fixed, impermeable type 
(the conventional use). 



Proof: Notice that an expression involving the kinetic energy can be formulated by 
substituting v for w in the variational momentum equation (|^). Then 

p(v,^j) = (pVv,F~'j) ^ . -2fi(D(v),b(i)j) ^ . 

\ ^ ^ / L2{n) 

+p(v,bj) +(v,PN) . (12) 

The order of integration and differentiation are fully interchangeable (the volume is still 
a material volume overall for the type of time-dependent limits associated with free 
boundaries). The term containing the pressure, that is 
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vanishes as a result of incompressibility (equation (^). Equation ([T^ ) can accordingly 
be rewritten 



d 



V J2 



rimposed 



L2(n) 



+p(v,bj) , - + (v,PN\ , ^ 
where Kf.. ^ , is zero for the present by virtue of the fully interchangeable orders of 

i imposed v i ^ ^ o 

integration and differentiation (its meaning will be made clear in the pages to follow). 
Using Lemmas |l| and ^ an expression 



dk{v) 
dt 



L2(r) 



rimposed V 



(13) 



is obtained, where K = -p 



Vj2 



L2{Q) 



is the total kinetic energy. 



The term F N- [v — v'^^^) vanishes at fixed impermeable boundaries since both F N ■ v 
and v"^^^ vanish under such circumstances (assuming the description becomes purely Eu- 
lerian there). This self-same term also vanishes at free boundaries according to Lemma 
^. Boundaries of an imposed velocity type need not be accounted for as a consequence 
of the stated "no forcing" condition, and so 



^ < -2uCK{v)+p{v:hj) 
dt - y J y\ ^ /^2(f^) 



(v,PN 



L2(f) 



This equation has a solution of the form 

K < k{v 

in the absence of forcing ( "no forcing" b ■ 



-2vCt 



PN = 0) 



A nonlinear, exponential-type energy dissipation in the absence of forcing is therefore an 
intrinsic property of the completely general referential description. This contractive flow 
property is also an intrinsic property of the conventional Navier-Stokes equations. 



3.2 Long— Term Stability under Conditions of Time— Dependent 
Loading 

The formulation of suitable load and free surface bounds is necessary before the issue 
of long-term stability (L^-stability) under conditions of time-dependent loading can be 
resolved. The energy transfer across boundaries at which there is an imposed velocity 
is a further factor which must be taken into account under conditions of forcing. The 
following lemma facilitates the formulation of load and free surface bounds. 
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Lemma 4 (Force, Free Surface Bounds) The inequality 

uC 



p (v, bJ 



L2(f2) 



v.PN) , < 



P 



L2(n) 



1 



2uC 



P 



L2(f2) 



+ 



L2{r) 
PN 



2 

L2(f) 



holds where vC is a constant, vC > 0. (Proof in Appendix IL) 



The relation immediately below will negate any convection-related contribution to the 
energy bound at imposed velocity-type boundaries. 



Lemma 5 (Convective Energy Rate at an Imposed Velocity-Type Boundary) 
The relation 

holds for boundaries at which there is an imposed velocity provided there is no nett in- 
flow/outflow across such boundaries and the description there becomes pure Eulerian. 
(Proof in Appendix IL) 



This done, the mathematical machinery necessary to the long-term stability analysis is 
in place. 



Theorem 2 (Long-Term Stability) A sufficient condition for the completely gen- 
eral reference description to inherit the property of long-term stability 

hm sup K{v) < —— 

under conditions of time-dependent loading, where this time-dependent loading, the speed 
of the surface and any imposed boundary velocity is bounded in such a way that 



P 



bj2 



PN 



L2(f) 



ij yi- ) i imposed V — 



is that the description is of a vanishing ii ■ {v — v'"^^^ type at free boundaries, that it 
becomes purely Eulerian at boundaries across which there is an imposed velocity or where 
boundaries are of a fixed, impermeable type. 



Proof: Cognizance must now be taken of a previously unencountered boundary type; 
that of a stationary boundary across which there is an imposed velocity. The limits of 
the integral on the left hand side of equation (|12D are time-dependant under such cir- 
cumstances and the volume is no longer a material volume overall. Kf,. ^ . in equation 

C3 i imposed V J- 

^The additional terms Kf,.^^^^^^^ - are given in Appendix I. They are only applicable in instances 
where there is an imposed velocity at the boundary. The other two boundary terms are only applicable 
at boundaries which involve tractions. 
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(p!3|) is no longer zero. Using Lemmas 0, H and || in equation ([131) , then applying the 
above bound, 

dk{v) ^ M2 

which, when solved, yields 

M2 



This in turn implies 

lim sup /^(i;) < 



2//2C2 ' 



The preceding analyses lead to natural notions of nonlinear dissipation in the absence 
of forcing and long-term stability under conditions of time-dependent loading for the 
analytic problem. These properties are also intrinsic features of real flows and the Navier- 
Stokes equations. 



The Energetic Implications of the Time Discreti- 
sation 



This section is concerned with establishing a class of time discretisations which inherit 
the self-same energetic properties (nonlinear dissipation in the absence of forcing and 
long-term stability under conditions of time dependent loading) as the analytic problem, 
irrespective of the time increment employed. In this section a generalised, Euler difference 
time-stepping scheme for the completely general reference equation is formulated and 
the energetic implications are investigated in a similar manner to that carried out for the 
analytic equations in the previous section. 

This stability analysis is inspired by the approach of others to schemes for the conventional 
Navier-Stokes equations. The desirability of the attributes identified as key energetic 
properties is recognised and they have been used as a benchmark in the analysis of 
various of the conventional, Eulerian Navier-Stokes schemes by a host of authors. Related 
work on the conventional, Eulerian Navier-Stokes equations can be found in a variety of 



references, for example Temam jl^ and SiMO and Armero [|T4 



The analyses presented here are extended, not only in the sense that they deal with 
the completely general reference equation, but also in that non-zero boundaries, so- 
called free boundaries and time-dependent loads are able to be taken into account (the 
former two as a consequence of the new inequality). The findings of this work have 
profound consequences for the implementation of the deforming reference equations. It 
is significant that many algorithms used for long-term simulation do not automatically 
inherit the fundamental qualitative features of the dynamics. 
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A Generalised Time— Stepping Scheme 



An expression for a generalised Euler difference time-stepping scheme can be formulated 
by introducing an "intermediate" velocity 

Vn+a = a V |t+At +(1 - a) V \t for a E [0, 1] (14) 

to the variational momentum equation (equation (^) on page |^ where v \t and v \t+/\t 
are the solutions at times t and t + At respectively, At being the time step. It is in this 
way that a generalised time-discrete approximation of the momentum equation, 

{pVw, KlJn-,a) ^^^^^^^^ - 2/. {dm, r>(i;„+.) J„+.)^,^-^^^^ 

+p(w,bn+aJn+a).^,^ + (w , P n+aN n+o) . ^.^ (15) 

is derived, where { ■ ) l2(^q ^ ^ denotes the inner product over the deforming domain 

at time t + a At. r„+a, F^+a, Jn+a, D^+a, Pn+a, and bn+a are likewise defined to be 
the relevant quantities evaluated at time t + aAt. 

It will presently become apparent that relevant energy terms are not readily recovered 
from the time-discrete equations for deforming references in general. It may therefore 
make sense to perform the analyses for the time-discrete equation in the context of diver- 
gence free rates of reference deformation only. A practically less restrictive alternative is 
too labour intensive. This investigation is accordingly restricted to a subclass of reference 
deformations in which "reference volume" is conserved. This is for reasons of expedience 
alone and it is hoped that this subclass of deformations is thought to be representative. 



Assumption 1 The assumptions Jn = Jn+a o-nd Jn+i = Jn+a O'l"^ made so that the 
desired energy terms are readily recovered as 



and 



K{Vn+l) = 



~i 


2 1 


_ 1 


2 


VnJn 


- = oP 


'-'n'Jn+a 













"^n+lJ-n+l 



'^n+lJn+a 



J J d J 

Remark: Notice that "^^^ — - = Jn+adivVn+a, the discrete form of — = Jdivt;'"'^-^, 

can consequently be rewritten as 

div vl^la = 

under the conditions of the above assumption. It is for the practical expedience afforded 
by Assumption |l| alone that this analysis is limited to instances in which divvl^l^ = 0. 

The following lemma will facilitate the elimination of the rate of energy change associated 
with the convective term under these conditions. 
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Lemma 6 (Discrete Convective Energy Rate) The following relation involving 
the discrete convective term holds for an incompressible fluid under circumstances of 



~P (vn+a, (VVri+a)-^^ n+a ('^n+a — "^n+a) "^"^ 

\ ^ ^ /L2(n„+„) 

= - 2^ ("^n+a, Vn+a {fJ_^^N ■ - V^+a)) J, 

(Proof in Appendix II.) 



Remark: Recall that in the investigation of the analytic problem, a term arising from the 
manipulation of the acceleration containing term (the term containing the rate of change 
of the Jacobian) cancelled with the convective energy. It is therefore not surprising that 
assumptions pertaining to the acceleration containing term (in particular to the rate of 
change of the Jacobian) in the discrete problem will, once made, also be necessary for 
the corresponding discrete convective energy term to vanish (reffering to the div f*"^-^ = 
condition of Lemma P). This is a good prognosis for the energetic behaviour of the 
discrete problem in circumstances of reference deformations excluded by Assumption 
The full ramifications of Assumption |l] are considered in Subsection |8.2| of Appendix I. 
This concludes the preliminaries required for the analysis of the time-discrete equation. 



4.1 Nonlinear Dissipation in the Absence of Forcing 



The following analysis establishes a class of time-stepping schemes which exhibit nonlin- 
ear dissipation in the absence of forcing regardless of the time increment employed. 



Theorem 3 (Nonlinear Dissipation in the Absence of Forcing) Suppose that 
the description is of a vanishing rin+a ■ {^n+a — ^n+a) type at free boundaries, purely 
Eulerian at boundaries across which there is an imposed velocity and that the deformation 
rate of the reference is divergence free. A sufficient condition for the kinetic energy 
associated with the generalised class of time-stepping schemes to decay nonlinearly 

K{vn+i) - K{vn) < -At2fi 

in the absence of forcing and irrespective of the time increment employed, is that the 
scheme is as, or more, implicit than central difference. That is 

1 



DiVri-Un) Jn 



'n+a)^n+a 



Proof: Expressing the "intermediate" velocities and Vn+a in terms of equation 



(|1^) and subtracting, the result 



Vn+a = {^-\) (^"+1 - ""n) + (16) 
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is obtained. The first step towards formulating an expression involving tlie kinetic energy 
of the generahsed time stepping-scheme (^) is to replace the arbitrary vector, w, with 
Vn+a- By further substituting (|T6|) into (|T5|) and eliminating the pressure containing term 
on the basis of incompressibility (equation (|^)), an expression involving the difference in 
kinetic energy over the duration of a single time step is obtained. 

Incompressibility and a restriction on reference deformations to those for which diwl^_l^^ 
is zero ensure that the Lemma |6| condition is satisfied. 



The equation 



-At 2fl 



~i 

2 



( "^n+a 1 P n+oi-^ n+a i -.^ ,:z. 



L2(f„+,) 



(17) 



is then obtained. The term F Jj^^N n+a 



{vn+a — Vnla) vauishcs at fixed impermeable 
boundaries since both F^"_^^Nn+a ■ Vn+a and v^^^ vanish under such circumstances (as- 
suming the description becomes purely Eulerian there). This self-same term also vanishes 
at free boundaries according to Lemma ^. Boundaries of an imposed velocity type need 
not be accounted for as a consequence of the stated "no forcing" condition, and so 



K{Vn+l) - k{Vn) < 



-At 2fi 



because of this condition. Thus the kinetic energy inherent to the algorithmic flow de- 
creases nonlinearly in the absence of forcing, irrespective of the time increment employed 
and for arbitrary initial conditions provided that 



a > 



and 



div -y^ta 



0. 



The former requirement translates directly into one specifying the use of schemes as, or 
more, implicit than central difference. Only for descriptions which are divergence free 
has it here been guaranteed that energy will not be artificially introduced by way of the 
reference. 

Remark: Notice (by Lemma P that for a = ^ an identical rate of energy decay 

k{Vn+l) - k{Vr 



At 



< 



is obtained for the discrete approximation as was obtained for the equations. 
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4.2 Long— Term Stability under Conditions of Time— Dependent 
Loading 

This second part of the time-discrete analysis estabhshes a class of time stepping schemes 
which exhibit long-term stability under conditions of time dependent loading, irrespective 
of the time increment employed. The following lemma is necessary to the analysis and 
is concerned with devising a bound for the energy at an intermediate point in terms of 
energy values at either end of the time step. 

Lemma 7 (Intermediate Point Energy) The following hound applies 

K{vn+a) >a{a-c + ac) k{vn+i) + {1 - a) (^1 - a - k{vn) 
where c is some constant, c > 0. 

The optimal choice of the constant c is established farther on. The following theorem 
establishes a class of time-stepping schemes which exhibit long-term stability under 
conditions of time-dependent loading regardless of the time increment employed. 



Theorem 4 (Long-Term Stability) Suppose that the description is of a vanishing 



n 



n+a 



i35 ^ 



n+a '-'n+a J 



type at free boundaries, that it becomes purely Eulerian at fixed 
impermeable boundaries and that the rate at which the reference is deformed is divergence 
free. A sufficient condition for the algorithmic flow to exhibit long-term stability under 
conditions of time-dependent loading assuming this time-dependent loading and the speed 
of the free surface is bounded in such a way that 





_1 


2 




p 


h 

"n+a'-' n+a 


+ 


P IV 

-'■ n+a-'- ' n+a 











+ V^C^\\V 



I 



IS 



1 

">2- 



Proof: Substituting Lemmas |], ^ and || into equation (|T^, applying the above bound 
and choosing a > | one obtains 

VL^K{Vri+a) < 



At ' ^-"^-"-"^ - 2uC 



From this point on the argument used is identical to that of SiMO and Armero [|14| for 
the conventional, Eulerian Navier-Stokes equations. Substitution of Lemma |^ leads to a 
recurrence relation, 

K{Vr^+l) < — -— ■ -j- K{Vn) 



1 + uCa{a- c + ac)At ' ' 2uC [1 + uCa{a - 1 + ac)At]' 



^Note that this bound does not incorporate a contribution from boundaries of an imposed velocity 
type in any obvious way. The two boundary terms are only applicable at boundaries which involve 
tractions. 
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Using this recurrence relation to take cognisance of the energy over all time steps, 
'1 - iyC{l - a)(l -a - ^)At' 



K{Vn 



+1) 



< 



1 + z/Ca(a — c + ac)At 
M^At 



Kivo) 



n-l 

E 



2z/C [1 + uCa{a - c + ac) At] [ 1 + i^Ca{a - c + ac)At 



[1 - uC{l - a){l - a - ^)At) 



(18) 



is obtained. An infinite geometric series which converges so that 

A/PAt 



lim sup K{Vn^ 



< 



1 



2uC [1 + uCa{a - c + ac) At] 

(l-i/C(l-a)(l-a-f)At) 
1 + i'Ca{a — c + ac)At 

M2 



2iyC [vCa{a -c + ac) + vC{l - a)(l - a - f ) 

results, providing the absolute ratio of the series is less than unity. That is 

l-i/C(l-a)(l-a-f)At 



1 + vCai^a — c + ac)At 



< 1. 



Therefore either 



a 



or 



-1 - ijCa{a - c + ac)At < 1 - ijC{1 - a) ( 1 - a ) At 



1 - iyC{l - a) [1 - a - - ) At < 1 + uCa{a -c + ac)At 



(19) 



in order for the bound to exist. Notice, furthermore, that for this desired convergence to 
be unconditional (regardless of the time increment employed) requires 

a-c + ac>0. (20) 

The denominator in the series ratio might otherwise vanish for some value of At. 
equation ( [T9| ) and equation ( pOD together imply 

(I -a) 



For a G 



< c < 



a 



a 



which in its turn implies 



:i-a) 



< 



[l-a] 



a 



a 



[l-aY 
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The choice of the parameter a > | therefore leads to an infinite geometric series which 
forms the desired upper bound. The minimum value of this bound occurs for c chosen 
according to 

1 1 



inf 



(i::ii)<c< " vCaia -c + ac) + vCil - aMl - a - ^) uC(2a - IV 

The value of this upper bound, which occurs for the choice of the parameter a > |, is 
then 

hm sup K('U„+i) < 

In this way one arrives at a class of algorithms which are unconditionally (irrespective of 
the time increment employed) stable. 

Remark: Notice that for a = 1 one obtains an identical energy bound for the discrete 
approximation as was obtained for the equations. 



5 Some Numerical Examples 



Some numerical results for problems of the type in question are presently given. The 
theory thus far developed was employed in the simulation of a driven cavity flow, a 
driven cavity flow with various, included rigid bodies, a die-swell problem and a Stokes, 
second order wave. 

The approach taken when approximating free surfaces, was that they may be treated 
as a material entity, that is, the material derivative of the free surface was assumed 
zero. Euler's equations and conservation of linear momentum were used to determine the 
motion of the rigid body. A predictor-corrector method was used to solve the combined 
sub-problems. 

A backward difference scheme was used to approximate the time derivative in the fluid 
sub-problem (in compliance with the Theorem |^ and Theorem ^ conditions) , the finite 
element method was used for the spatial (referential "space") discretisation and a 
element pair was used as a basis. A penalty method was employed to eliminate pressure 
as a variable and nonlinearity was circumvented by way of a new, second order accurate 
linearisation. Linearising with a guess obtained by extrapolating through solutions from 
the previous two time steps leads to second order accuracy. 



Theorem 5 The linearised terms, {2v \t —v \t-At) ■ Vf \t+At and v \t+At ■V(2t; \t 
—V \t-At), second order accurate (have error O(At^) ) approximations of the nonlinear 
term {v ■ Vt») \t+At- 



Proof: 



V \t+At 



V 



. dv 



0{At^ 



(by Taylor series) 
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[v ■ Vv) \t+At 



V I, +At I "^'^ ^J*-^* + 0(At) j + OiAt') 



'2v \t -V \t-M +0{At^) 
2v \t -V +0{At^] 



(using a backward 
difference) 



Vv 



t+At 



[2v \t -V \t-At] ■ \t+At +OiAt^ 



The above linearisation schemes are an improvement on the conventional v \t -Vv \t+At 3 
or V \t+At 'ViJ \t linearisation schemes by an order of magnitude. A detailed exposition 

of all numerical methods otherwise used in these simulations can be found in Childs 

and Reddy HI. 



5.1 Example 1: Driven Cavity Flow 

The problem is essentially that of a square, two-dimensional pot whose lid is moved 
across the top at a rate equal to its diameter for a Reynolds number of unity. The 
boundary conditions are accordingly "no slip" on container walls and a horizontal flow 
of unity across the top (depicted in Fig. |]). 




6 The Leaky Lid Problem 1 The uniform mesh (144 elements). Navier-Stokes pressure (control). 



Figure 3: The Problem, the Mesh and the Pressures Obtained Using the Conventional 
Eulerian Equations. 



The idea here was to compare results obtained using the completely general reference 
equation on a deforming mesh with those obtained using the conventional, Eulerian 
Navier-Stokes equations. 

^Favoured in terms of both rate and radius of convergence by Cuvelier, Segal and van Steen- 

HOVEN §. 
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Pjmjje fbi a 5% ITS Ji compioE- The d'ElalediraJi (I'Uelemen'fc] PjoEUje fbj a 5% maJi decompue- 
EEiDn. "competed" by 3%. mion. 



Figure 4: Pressures Obtained Using the Completely General Reference Equation and a 
Deforming Mesh. 



The corresponding velocity profiles along the cuts depicted in Fig. |^ are given in Fig. ^. 
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Velocity Profiles 




Figure 5: In this test part of the mesh was successively compressed and decompressed by 
5 % over two time steps of length 0.05. 
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5.2 Example 2: "Pebble in a Pothole" 

In this example rigid bodies of varying mass and moments of inertia were released from 
rest in a flow dictated by the same boundary conditions as the driven cavity flow of the 
previous example. One would expect a die bead (a small rigid body of neutral bouyancy) 
to move in tandem with the fluid soon after its release from rest. One might also expect 
a clockwise rotation to be induced by concentrating the mass closer to the centre i.e. 
lowering the moment of inertia. 

The finite element mesh was automatically generated and adjusted about the included 
rigid body in what is possibly a slightly novel fashion. A small region of mesh immediately 
adjacent to the included rigid body was repeatedly remapped to cope with the changing 
orientation, the remainder was squashed/stretched according to the translation. 

To begin with, a square region of mesh centered on, and including the rigid body, is 
deleted (depicted in Fig. |^). Each of four wedge-shaped regions is then demarcated (the 
intersections of lines which bisect corners and edges of the square frame, with the surface 
of the rigid body are located using Newton's method) by as many points as there are 
nodes in an element i.e. each wedge shaped-region is set up as a massive element. 




(-1. 1) (1. 1) 




(-1. -1) (1. -1) 

Figure 6: The Local Distortion is Obtained by Mapping Square Chunks of Rectangular 
Mesh Using Finite Element Mappings. 

Chunks of uniform mesh, which have identical extremities to those of the master element, 
are then mapped into the newly-demarcated, wedge-shaped regions using finite element 
mappings (in exactly the same manner as points in the master element domain are. 
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in theory, mapped into individual mesh elements). Further, fine adjustment of nodes 
intended to delineate the surface of the rigid body is accomplished by moving them 
along a line between node and centre, to the rigid body surface using Newton's method. 
The mesh outside the "box" (the box containing the 4 wedges enclosing the rigid body) 
is squashed/stretched according to the requirements of the translation (the nodes are 
translated by a factor inversely proportional to their distance from the box). This method 
satisfies the requirement that h ■ (v — v'^'^^^ vanishes at the fluid-rigid body interface (a 
condition in Theorems |I], ^ |^ and ^). Mesh refinement in the vicinity of the included, 
rigid body is an automatic by-product of this method. 





Figure 7: Typical meshes which result when using this method of automatic mesh gen- 
eration about rigid bodies which are simultaneously rotating and translating. 



Various rigid bodies were introduced to the driven cavity fiow problem described in 
Subsection Ol in the absence of a body force. The results in Fig. involve the ellipse 



^2 ™2 

I + ft = 0.025^ 



whose major axis is 0.1. The quantities m and JM{no sum) are a dimensionless mass and 
ith principle moment of inertia respectively. 
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Figure 8: The trajectories of various included rigid bodies released from rest at the centre 
of the driven cavity flow described. Top Left: Re = 0.025, rh = 251.3, J33 = 314.2 and 
t = 3.6 sees. Top Right: Re = 0.025, rh = 251.3, J33 = 1.0 and t = 4.0 sees. Bottom 
Left: Re = 0.025, m = 251.3, J33 = 0.1 and t = 3.6 sees. Bottom Right: Re = 1, 
rh = 1, moment of inertia (scaled) = 0.1 and t = 2.0 sees. 



5.3 Example 3: Die Swell Problems 

The axis-symmetric die swell (or fluid jet) problem is a free surface problem well doc- 
umented in the literature (Kruyt, Cuvelier, Segal and VAN Der Zanden |||], 
Omodei [|12| and Engelman and Dupret quoted in Kruyt et al.). The basic theme 
to this problem is the extrusion of a fluid with initial parabolic flow profile from the end 
of a short nozzle. 
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Figure 9: Die swell ratios predicted for various Reynolds numbers using an inlet velocity 
profile of Vi = — x|) and the methods described. (Bars on the variables merely 

indicate that they are dimensionless.) 
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Figure 10: Die swell ratios predicted for various Reynolds numbers using an inlet velocity 
profile of vi — ^^{1 — xl) and the methods described. 
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Figure 11: Die swell ratios predicted for various Reynolds numbers using an inlet velocity 
profile of — ^ ' 



2Q "~ ^2) ^^"^ methods described. 



5.4 Example 4: A Stokes Second Order Wave 

In this problem the velocity profile and surface elevation predicted by Stokes second 
order wave theory (see Koutitas [^) were used as boundary conditions for flow and free 
surface subproblems respectively. 




Figure 12: A Stokes second order Wave. 
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Problems with wave propagation were subsequently experienced as time progressed. It 
should be noted, however, that the problem was not attempted with the same seriousness 
as previous examples and the mesh was poor (there being only three elements in the 
vertical extent of the mesh). 



6 Conclusions 



The correct equations, which describe the motion of an incompressible, Newtonian fluid 
and which are valid for a completely general range of reference deformations, are equations 
(|) and (|^. For implementations requiring the equations to be evaluated about a single 
instant within each time step only (eg. finite differences), the deformation gradients may 
be assumed identity i.e. the equations of Hughes, Liu and Zimmerman [^] (equations 
(0) and (|11])) will suffice. 

In this work it is shown (as was hoped) that nonlinear, exponential-type dissipation in the 
absence of forcing and long-term stability under conditions of time dependent loading are 
properties automatically inherited by such deforming reference descriptions. The single 
provisor is that the conventional boundary descriptions are used (vanishing h- (^v — v'^'^-^^ 
type at free boundaries, purely Eulerian at boundaries across which there is an imposed 
velocity or where boundaries are of a fixed, impermeable type). These properties are 
intrinsic to real flows and the conventional, Eulerian Navier-Stokes equations. 

Relevant energy terms are, however, not readily recovered from the time-discrete equa- 
tions for deforming references in general. Only for divergence free rates of reference 
deformation could it consequently be guaranteed that energy would not be artificially 
introduced to the algorithmic flow by way of the reference. A further casualty of the dis- 
crete analysis is its failure to account for flows driven by their boundaries in any obvious 
way i.e. boundaries of an imposed velocity type do not enter explicitly into the bound. 
Scope for the further development of this work therefore exists. 

The divergence free assumption was made for reasons of expedience alone and it is hoped 
that the findings of the time-discrete analysis can be extrapolated to a more general 
class of mesh deformations. If one were to be overly cautious on this basis one would be 
faced with the additional challenge of enforcing mesh deformations which are divergence 
free. Such a totally divergence free description may, however, not be practical. Both the 
purely Lagrangian and purely Eulerian fluid descriptions have divergence free rates of 
distortion. 

What is clear is that there are inherent problems with using certain classes of time- 
stepping schemes and the use of finite difference schemes more implicit than central 
difference is consequently advocated. The limitations of the time-discrete analysis do not 
detract from this finding in any way. Such differences exhibit the key energetic properties 
(nonlinear, exponential-type dissipation in the absence of forcing and long-term stability 
under conditions of time dependent loading) irrespective of the time increment employed. 
A backward difference is the obvious choice. Calculations at time t+aAt would require an 
intermediate mesh and associated quantities for instances in which a 7^ 1 (since a > 
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The author recommends a strategy in which a predominantly Eulerian description is 
used, where possible, for the bulk of the problem (from an efficiency point of view) and 
the completely general reference description for the remainder is appropriate. Purely 
Eulerian descriptions have the advantage of a "one off" finite element construction and 
involve none of the hazards of a badly distorted reference. 

With regard to numerical implementation and using a element pair, it was found 

that pressures approximated as linear on the master element still led to so-called "lock- 
ing" or "chequerboard" modes. The pressures needed to be linear on the actual elements 
themselves. This finding makes sense if one considers that a linear function mapped from 
the master element using a Q2 mapping will no longer be Pi for non-rectangular ele- 
ments (the element pair was shown to satisfy the L.B.B. condition in the context 
of rectangular elements). 

Lastly, the linearised terms {2v \t —v \t-At) ■ Vt; \t+At and v \t+At ■V(2i; |( — 1> \t-At) are 
second order accurate approximations of the convective term and a remarkably practical, 
simple and effective method to automatically generate meshes about included rigid bodies 
was devised. 
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8 Appendix I 



8.1 The Contribution K 



rimposed 



at Imposed Velocity— Type Bound- 



aries 



Additional terms which arise from the limits of the integral 



P 



T.0 lz.O Z.0 
1 -^2 -^3 



at 



when changing the order of differentiation and integration at boundaries across which 
there is an imposed velocity are: 



1 



^0 J^o 



V ■ V J (1x^(1x2 



Dxl 

15t 



V ■ v,Jdx^dx2 
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+ 



+ 



*i Dx'o 



Dt 



V ■ vJdx^i 



dxi 



x? Dt 



V ■ vJdx?, 



dxi 



xO JxO Dt 



V ■ vJ 



dxodxx — / — (v ■ vj] „ dxodxi 
's JxO JxO Dt ^ ) xl ^ ^ 



ix\ JxO Dt 



(21) 



(using Leibnitz's rule for differentiation under the integral sign). 
Note that the terms „ . can be given a 



Dt 



d... 



^-1 , 



+ V . . . F {v-v 



interpretation in terms of the relation (|l]) established at the beginning of Section |2.2| . For 
descriptions which are purely Eulerian at such fixed boundaries, v'^'^^ vanishes and F is 
identity. Thus the extra terms, (|2T| ) above, can be rewritten (with minus sign omitted) 



2^ 



1x0 ^ 



"-2 -^3 

4 



Vx'i ■ V i I V ■ V dx'idxo I — Vx? ■ v i j v ■ v dx^dxo 



V ■ V dx^ 







dxi — 



/•x[ / rx'^ 

JxO ^''''''[jxO ''■''^'^^^ 



dxi 



+ / Vx'n ■ V (v ■ v)\-, dx2dxi - / / Vx^ ■ V (v ■ v)\~o dx2dxi 

JxO JxO '""s J^o J^o 

i.e. a total rate of energy transport across the boundaries, similar, but not identical to 
/p(n ■ v){v ■ v)dV. These terms, dubbed -^^fimposod i work, must be added to the 

right hand side of equation (|13|) when circumstances require. 

This approach may seem comparatively crude in the light of rather elegant work done by 
Temam [0 for flows driven by their boundaries, however, the intended purpose differs 
slightly. What is here being sought is a bound formulated in terms of known, physically 
comprehensible quantities at the boundary which are independant of the solution. 



8.2 The Ramifications of Assumption |I] 



If one does not make Assumption |I], contributions from the 



term in equation ( |15D and 



1 / ~ ~ {Jn+l " Jn) 



At 



in Lemma |^ amount to 
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P 1 



At 


[a 




1) 


P 




1 




At 


[a 




I) 


P 




1 




At 


[a 




\) 



'^n+oi '^n+i '^'i+a/ 2 



3a^— (K{v^+,) - k{v. 



At 



- {oLVn+i + (1 - a) w„, aVn+1 + (1 - a) VnaJn+i + (1 - a) J„ 
- (avn+i + (1 - a) ^ + Vn) aJn+i + (1 - a) J„ 

- ^ + (1 - a) -^n, Vn+a^ + J„ 



+ 



2 - 3a) {vn+i, Vn+iJn) + a(2 - 3a) i;„+i J„,+i) 



At 

1 / _ (1 — a)^(6a — 1) 
--(1 - a) (3a - 1) v„+iJ„,) H _ J„+i) 



(by repeated substitution of ([l^ ) and (|T6|)). Thus the ramifications of Assumption |l| are 
that the total 



-^6 (a--\kiv^ + 4- 
At \ 2) ^ ' At 



a 



(2 - 3a) {Vn+l, Vn+lJn) + "(2 - 3a) Vn+lJn+l) 



\/o IN /- ~ T\ (1 — a)^(6a — 1) „ ^ , 

-|(1 - a)(3a - 1) {Vn,Vn+lJn) H _ {Vn,VnJn+l) 



is positive, furthermore it is sufficiently positive to offset any subsequent short-coming 
which arises when Assumption |1| is exploited in the proof of Lemma |^ i.e. in the event of 



and 



not being positive. 



9 Appendix II (Proofs) 

Proof of Relation n 



The above relation (taken from Hughes, Liu and Zimmerman 0) is obtained by re- 
calling that the material derivative (total derivative) is the derivative with respect to 
time in the material configuration. Thus 
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dt dxj dt 



(22) 



A more practical expression is needed for (the velocity as perceived in the distorting 
reference). This can be obtained by considering 

Xk{xo,t) = Xl{X{xQ,t),t) 



so that 



dX, 



dt 



dXl 



or 



X(i fixed 



dxj I dXk 



dt 



+ 



see Fig. |1]) 

dXl dX, 



X fixed 



ax, 



dt dxk \ dt 



dXl 



dt 



X fixed J 



Xq fixed 

Substituting this expression into equation (p2D, the desired, suitably practicable result is 
obtained. 



Proof of Inequality |1] 

Consider the change to spherical coordinates 

'Uj(r, 6*, (p) = Vi{r sin 6 cos <!) — x™'^", r sin 6 sin — x™'^", r cos 6 — x™'^") 

centred on a;°"^™. Suppose the radial limits of the domain and neighbourhood are denoted 
Rb{0, (f>) and Ra{0, (j)) respectively. By the fundamental theorem of integral calculus 

/ \ 2 / /■'" dv- \ ^ 



jRj9,<h) £ dr 



2 



lRaie,<f>) ^ dr 



2 



jRaie,<p) £^ JRa(e,<i>) \dr I 



(by Schwarz inequality) 



2 



/■Rmax 1 rRbid,4>) / dVi I o 

- Ir T^'^ Lie,, b'^^'^''^) ^^'^^^^^^ 

"'^min S JRa{e,,) \Or I 



\ J- •'max -"'minj / ^ ' / ^fJi 
(^Rniax -^min) 



jRa(e,d>) \ dr / 



Rmax Rmm 

2 



where \4(^,0)= / -^(^,^,0) e'^^ 

jRaie^ii) \dr I 
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Integrating this result over that part of Q outside the neighbourhood (angular extent 
being e„(0) <e< Qbicj)) and < < ^b) 

/ / [Vi{r,9,(f)) - Vi\R^(0^^)] r sin edrd9d(j) 



-Rmax-Rmin -^^o ''Qai't>) JRa{d,4>) 
(B —B) r^b r^b{4>) V / /-^Jmax „ \ 



3i?max^min J^a J@a{<t>) JRa{9,<f>) \ Or 



^ I, -"-ma x -"-mm A -'''max -"-min^ / / / 



3-Rmax-Rmin -^^a •^0o(</') JRa(d,4>) 



1 / \ ^ 



~'~r^sin^ 6 [d(l) 



dr I r"^ \d9 J 
sin 9drd6d(p 



_ p . V7?3 _ 7?3 /•e6(<^') /•Ki.ce, 



/•y^t r'^b[(f') fnb(<^<<P) 

/ / / (Vvi) ■ (Vvi)r^ sin edrdedcf). 

J^a Jea(<t>) JRa(e,4>) 



3 Rmax -^min 

Changing back to the original rectangular coordinates and defining v |bndry to be a radially 
constant function throughout Q which takes the values of v for r = Raid, 0), 

/ (Viix) - V, \^^,^f dn < (-^'na.-^min)(i^Lx-i^iJ / ^^y^^a^)) ■ {Vv,{x)) dil 

where is Q excluding the neighbourhood. Summing over i, 

I {V-V\^^^,,)-{V-V\^^^,,)dn < (^"^ax-^min)(gLa.-^^in) f (y^) ^ (y^) 

Making use of either the Cauchy-Schwarz or triangle inequality, 

-.ill 11-1,1 II ^^ ^ (-^max ~ Rmin){R^ax ~ -^min) ||V7^.||2 

^llL2(n,) IF Ibndry||i2(f^^)j ^ P 1 1 V 1 1 L2(n^) , 

^ ^max -'•'min 



and remembering that sup | v |ij„(e,0)|< c, 

(-^^ax ~ Rmhi) {Rmiix ~ RminJ 



I^IIl2(q*) - 



'^RmnxR] 



max-* i-min 



Consider the terms ||t^||^2 and ||c||^2. Comparing these terms under circumstances of 
sup I V I < c leads to the conclusion that the inequality holds over the neighbourhood and 
that the inequality is therefore unaffected when the domain of integration is extended 
to include the neighbourhood. Of course, the radial extension of v |bndry can be used in 
place of c in instances where inclusion of the neighbourhood is not required. 
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Proof of Lemma |l] 



If, in particular, v 



bndry 



in Inequality |l|. 



L2(n) 



C 



ll^llL2(f7) < \\D{v)\\l,^^^ 

(The relationship between D and arises in the context of the original equations 
involving diver. It is because 



D 



2 ("^^j'i "I" "^ijO 



2'^«Ji 



(changing the order of differentiation) 
(clivt? = by incompressibility) , 



assuming, of course, that v is continuous and differentiable to first order.) Rewriting in 
terms of Cl 



-kiv) < 
p 



D(v)j2 



Proof of Lemma |2| 



^-1 



Consider u ■ {'Vv)F wJ: 

UiVijFji}wkJ = -UijViFfj^WkJ - UiVi{Fjj}wkJ)j + {uiViFjf}wkJ),j 

by the product rule. In the terms arising from {F~f}wkJ)j, both F~f}j and J^Ff^} van- 
ish under the condtions specified (in section P^ ) for equations (|TD|) and ( [TT| ) to be a 
completely general reference description. Thus 

UiVi^jFff}wkJ = -UijViFf^^WkJ - UiViF-f}wkjJ + {uiViFfj}wkJ),j. 

Integrating over the domain Q and applying the divergence theorem, 

V, {\/u)F'^wJ' 



u,{\/v)F ^wj) 



L2(n) 



u[Vw:F ^] ,vJ 



+ {u,v{F-'n.w)j)^,^^^^ 



Thus the term 



2{v,{Vv)F ^ (v-v^'f) j) 



I L2(fi) 



v{V [v- v^"^) : F 



vJ 



v,v(F ' N ■ {v - J 



L2(f) 



(23) 
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v(Vv''^:F ),vJ, 



L2(f) 



(by incompressibility) 



v,v(F ' N ■ {v - i'^'f)) J 



(v 



L2(f) 



dJ 



since — = Jdivv''^^ (which is JVv'^^^ : F ^) m the same vein as Jq = Jq div v (the 



kinematic result used in Section 



Equation (|23| ) is vital to the deforming reference analysis in particular. It forms the basis 
to this lemma and another (Lemma ^ concerned with the time discrete analysis. 



Proof of Lemma 



By referring to Fig. |I| one observes that 

\j[xQ,t) = X*{\{xo,t),t) 

and consequently that 



dt 



That is 



Xq fixed 



^-1 . 



dt 



+ 



X fixed 



dX* dK 
dxi dt 



dX 
'dt' 



— 1 _ ~ f 

In other words F {v — v^'^^) is the velocity perceived in the deforming reference. This 
perceived velocity is tangent to the free surface since the description was stipulated to be 
one in which h ■ (^v — v^^^^ vanishes at free surfaces. Remembering that iV is a surface 
normal as defined in terms of this self-same reference, 



(y — v'' 



at free boundaries. 
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Proof of Lemma 



In terms of the Cauchy-Schwarz inequality, 



v,bJ) - < 

L2(n) - 

< 



Vj2 

uC 



L2(n) 



+ 



L2(n) 



bj2 



for uC>0 



by Young's inequality. Similarly 



PN 



L2(f) 



for uC > 0. 



Proof of Lemma |5| 

For a description which becomes purely Eulerian at a fixed boundary across which there 
is an imposed velocity, 



~-p (^v, v{f 'n -{y- v''^)) J 



L2(f) 



1 



- -//^2(p) 



p{v,v{n-v))^2. 



< -p 



V ■ v{n ■ V 



(n ■ v)^ dr 



< -p { n ■ V dT 



2 ' ' / \2 ^2 

[v ■ V) n ■ V dV 



(by Schwarz inequality) 



_ ... flow .o.na.... X „ . Wr OU B.ea vo..e on_^^^ 



fiuid must be zero. 



Proof of Lemma ^ 

By equation (^) of the Lemma H proof. 



Proof of Lemma 



By Young's inequality 



n+a 



11 



< 'I 



+ 



2c 



2 



(24) 
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for c > 0. Writing K{vn+a) explicitly in terms of the "intermediate" velocity definition, 
(0), leads to 



k{Vn+a) = a'^k{Vn+l) + {1 - afK{Vn) +2a{l - a) (^Vn+l,VnJn 

> a^kivn+i) + (1 - ayk{v^) 
—2a(l — a 







_1 















> a [a — (1 — a)c] K{v„^i) + (1 — a) 



a 



[1-a)-- 
c 



K{Vn) 



using equation 
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